{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "try:\n",
    "    import qaths\n",
    "except ImportError:\n",
    "    # If we can't find qaths, then try to add some directory where it could be.\n",
    "    # First, we determine the current notebook directory.\n",
    "    import os\n",
    "    if not 'workbookDir' in globals():\n",
    "        workbookDir = os.getcwd()\n",
    "    \n",
    "    def is_root(path: str) -> bool:\n",
    "        path = os.path.realpath(path)\n",
    "        return path == os.path.dirname(path)\n",
    "    \n",
    "    # Then we try to find qaths in a parent directory.\n",
    "    import sys\n",
    "    current_dir = os.path.dirname(os.path.abspath(workbookDir))\n",
    "    import_successfull = False\n",
    "    # Loop in parent directories until we can import qaths or we \n",
    "    # find the root.\n",
    "    while not (import_successfull or is_root(current_dir)):\n",
    "        sys.path.append(current_dir)\n",
    "        try:\n",
    "            import qaths\n",
    "        except ImportError:\n",
    "            # Remove the added directory from the PYTHONPATH in order to\n",
    "            # not pollute it with a lot of useless directories.\n",
    "            sys.path.pop()\n",
    "            current_dir = os.path.dirname(current_dir)\n",
    "        else:\n",
    "            print(\"Found qaths library in {}.\".format(current_dir))\n",
    "            import_successfull = True\n",
    "    \n",
    "    # If qaths has not been successfully imported, warn the user.\n",
    "    if not import_successfull:\n",
    "        print(\"qaths not found!\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Imaginary number permutation matrix simulation\n",
    "\n",
    "This notebook illustrates how to use the `simulation.simulate_imaginary_float_weighted_hamiltonian` function to simulate $e^{-iHt}$ where $H$ is a permutation matrix with imaginary weights. In other words, $H$ must have at most one imaginary number in each row/column.\n",
    "\n",
    "### Examples of matrices $H$ concerned:\n",
    "\n",
    "A $\\cdot$ represents a value of $0$.\n",
    "\n",
    "$$\n",
    "H_1 = \\begin{pmatrix}\n",
    "\\cdot{} & i\\pi    & \\cdot{} & \\cdot{} \\\\\n",
    "-i\\pi   & \\cdot{} & \\cdot{} & \\cdot{} \\\\\n",
    "\\cdot{} & \\cdot{} & \\cdot{} & \\cdot{} \\\\\n",
    "\\cdot{} & \\cdot{} & \\cdot{} & \\cdot{} \\\\\n",
    "\\end{pmatrix}\n",
    "$$\n",
    "\n",
    "### Note on precision\n",
    "\n",
    "The imaginary numbers are currently represented as fixed-point numbers on qubits: $r$ qubits are dedicated to the integer part and $2r$ qubits are encoding the fractional part. The approximated value $\\tilde{H_{ij}}$ of an entry $H_{ij}$ is given by the integer value represented by the $3r$ qubits divided by $2^{2r}$.\n",
    "\n",
    "This representation introduce a truncation error of $\\mathcal{O}\\left( 2^{-r} \\right)$ in the resulting state."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Generate the program\n",
    "\n",
    "The following cell generates the quantum program that will simulate the given Hamiltonian. The variables `n`, `time`, `r` and `matrix` can be changed by the user to simulate different dynamics."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy\n",
    "import scipy.linalg\n",
    "import scipy.sparse\n",
    "\n",
    "from qat.lang.AQASM import Program\n",
    "\n",
    "from qaths.common.simulation.base.fixed_point_weighted import simulate_imaginary_float_weighted_hamiltonian\n",
    "from qaths.common.generation.fixed_point_weighted import generate_1or0_sparse_float_weighted_hamiltonian, matrix2float_weighted_permutation\n",
    "from qaths.common.qram import QRAM\n",
    "\n",
    "n = 1\n",
    "r = 3\n",
    "time = 1\n",
    "\n",
    "matrix = generate_1or0_sparse_float_weighted_hamiltonian(n, p=1, max_amplitude=2**r-1, imaginary=True).tocsc()\n",
    "matrix = scipy.sparse.coo_matrix(numpy.array([\n",
    "    [0, -0.681140613516],\n",
    "    [0.681140613516, 0]\n",
    "])).tocsc()\n",
    "permutation, weights, signs = matrix2float_weighted_permutation(matrix, r)\n",
    "matrix = 1.j * matrix\n",
    "\n",
    "M = QRAM(permutation)\n",
    "W = QRAM(weights)\n",
    "S = QRAM(signs)\n",
    "\n",
    "prog = Program()\n",
    "x = prog.qalloc(n)\n",
    "m = prog.qalloc(n)\n",
    "w = prog.qalloc(3*r)\n",
    "p = prog.qalloc(1)\n",
    "s = prog.qalloc(1)\n",
    "p2 = prog.qalloc(1)\n",
    "\n",
    "print(\"Applying... \", end=\"\")\n",
    "prog.apply(simulate_imaginary_float_weighted_hamiltonian(M, W, S, n, r, time), x, m, w, p, s, p2)\n",
    "print(\"Done!\\nConstructing circuit... \", end='')\n",
    "circ = prog.to_circ()\n",
    "print(\"Done!\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Visualisation is disabled by default because the circuit takes a very long time to render.\n",
    "# If the circuit is not too big, you can visualise it by uncommenting the line below.\n",
    "#%qatdisplay circ"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Simulate the quantum circuit and recover the amplitudes\n",
    "\n",
    "The following cell launch the linalg simulator to simulate the quantum circuit produced in the previous step.\n",
    "\n",
    "The MPS and Feynman simulators cannot be used for the moment as the QRAM gate may produce gates with an arity over 3.\n",
    "\n",
    "The statevector simulator is not used because it will store in memory a vector of $2^{3n+1}$ amplitudes instead of the $2^n$ amplitudes used by the current implementation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from qat.core.task import Task\n",
    "#Running on linalg.\n",
    "from qat.linalg import get_qpu_server\n",
    "\n",
    "task = Task(circ, get_qpu_server())\n",
    "\n",
    "task.execute()\n",
    "\n",
    "states = list(task.states())\n",
    "amplitudes = numpy.zeros((2**n, ), dtype=numpy.complex)\n",
    "for state in states:\n",
    "    if state.probability > 1e-10:\n",
    "        amplitudes[state.state._int_state] = state.amplitude"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Compare the results of the simulation with the analytical ones\n",
    "\n",
    "As we are still in a reduced dimension, classical computers can compute the exact result we expect to recover from Hamiltonian simulation. These results are computed with explicits tensor products and matrix exponentiation and then compared with the results obtained by simulating the quantum circuit.\n",
    "\n",
    "Warning: if one or more quantum gates are executed before the Hamiltonian simulation routine, `phi_0` needs to be changed accordingly. `phi_0` is the initial state of the simulation, `phi_t` is the final state."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "zero_state = numpy.array([1, 0])\n",
    "H_state = numpy.ones((2,)) / numpy.sqrt(2)\n",
    "phi_0 = 1\n",
    "for _ in range(n):\n",
    "    phi_0 = numpy.kron(phi_0, zero_state)\n",
    "exp = scipy.linalg.expm(-1.j * matrix * time)\n",
    "phi_t = numpy.dot(exp.todense(), phi_0)\n",
    "\n",
    "if numpy.allclose(phi_t, amplitudes):\n",
    "    print(\"Amplitudes are matching!\")\n",
    "else:\n",
    "    distance = numpy.linalg.norm(phi_t - amplitudes)\n",
    "    if distance < 1/2**r:\n",
    "        print(\"Amplitudes are matching with an error of {}.\\nAnalytic:\\n\\t{}\\nSimulated:\\n\\t{}\".format(distance, \n",
    "                                                                                                       phi_t, \n",
    "                                                                                                       amplitudes))\n",
    "    else:\n",
    "        print(\"Wrong amplitudes (2-norm distance of {}):\\nAnalytic:\\n\\t{}\\nSimulated:\\n\\t{}\".format(distance, \n",
    "                                                                                                    phi_t, \n",
    "                                                                                                    amplitudes))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.4.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
